/*=====================================================================================

 news_sni.do: 
	
   Creates the expected dividend based on APFC market data in APF_Data.xlsx	
   

		Lorenz Kueng, Feb 2015

=======================================================================================*/

cap log close PFD_03
log using "$homedir/log-files/PFD_03_$date.log", replace name(PFD_03)					


*================================================================================
* News Shock: Expected next dividend based on APFC data and public formula
*================================================================================

import excel "$homedir/data/APF_Data.xlsx", sheet("Monthly SNI") cellrange(A14:I310) clear 

	keep A B I
	rename A year
	rename B month
	rename I sni
	destring year,  replace
	destring month, replace

	generate Month = ym(year,month)
	format Month %tm
	order  Month
	tsset  Month

	drop if year<1991 // drop first few imputed observations


	* Deflate

	gen date = Month
	merge 1:1 date using "$homedir/data/stata/CPIu_monthly.dta", keepusing(cpiu)
	drop if _merge==2
	drop _merge date

	generate SNI_nominal = sni
	replace  sni = sni/cpiu*100

	tsset Month
	tsline sni SNI_nominal 


	*------------------------------------------------------
	* Unit-root tests.   Conclusion: no unit root -> d=0
	*------------------------------------------------------
/*
	* Graphically
	
	tsline sni

	
	* Augmented Dickey-Fuller
	*
	* -> reject null of unit root, hence stationary

	dfuller sni if year>=1991 & year<=2013, regress
	dfuller sni if year>=1991 & year<=2013, trend regress
	dfuller sni if year>=1991 & year<=2013, drift regress
	dfuller sni if year>=1991 & year<=2013, noconstant regress
	
	dfgls   sni if year>=1991 & year<=2013, notrend 
	dfgls   sni if year>=1991 & year<=2013, ers 
	
	
	* Phillips-Perron  ("Similar to Augmented Dickey-Fuller, but using Newey-West SEs")
	*
	* -> reject null of unit root, hence stationary

	pperron sni if year>=1991 & year<=2013, regress
	pperron sni if year>=1991 & year<=2013, trend regress
	pperron sni if year>=1991 & year<=2013, noconstant regress
*/


	*------------------------------------------------------
	* Model selection.   Conclusion: p=1, q=1
	*------------------------------------------------------

	* Non-parametrically (graphically) 
/*	
	ac  sni if year>=1991 & year<=2013, level(99)  // correlogram -> suggest q=1 or 6
	pac sni if year>=1991 & year<=2013, level(99) // partial correlogram -> suggests p=1 (or 5)
	* => arima(1,0,1) seems resonable
*/
	
	* ARIMAFIT: AIC, BIC, Log-Likelihood  -> suggest arima(1,0,1) or arima(8,1,8)
/*
	cap n net search arimafit
	cap n net from http://fmwww.bc.edu/RePEc/bocode/a/ /*Chris Baum's website*/
	cap n net install arimafit, from(http://fmwww.bc.edu/RePEc/bocode/a/)

	eststo clear
	forvalues p=0(1)10{
	 forvalues q=0(1)10{
	  display "p=`p', q=`q'"
	  eststo p`p'q`q': capture quietly arima sni if year>=1991 & year<=2013, arima(`p',0,`q') // fix estimation sample for future replications  
	  capture quietly: arimafit if year>=1991 & year<=2013
	  global AIC`p'`q'=r(aic)
	  global SIC`p'`q'=r(sic)
	  global LLF`p'`q'=r(llf)
	 }
	}
	est stat p*q*
*/
	
	
	* VARSOC: AIC, SBIC, HQIC, FPE, LR criterion  -> AR(5 or 6) is preferred
/*	
	varsoc sni if year>=1991 & year<=2013, maxlag(20)
	varsoc sni if year>=1991 & year<=2013, maxlag(20) noconstant
 	varsoc sni if year>=1991 & year<=2013, maxlag(20) lutstats
*/


	* Forecast Error Variance  -> most models perform similarly
/*
	cap drop time
	generate time = _n
	tsset time

	forvalues q=0(1)10 {
	forvalues p=0(1)10 {

		display _n(1) "arima(`p',0,`q')
		
		quietly{
		capture{ 

			forvalues i=1(1)12{
				cap drop sni`i'
				generate sni`i'=.
			}

			arima sni if year>=1991 & year<=2013, arima(`p',0,`q')

			sum time
			local T=`r(max)'
			forvalues t=1(1)`T'{
			
				cap drop temp
				predict  temp, y dynamic(`t')
				
				forvalues i=1(1)12{
					replace sni`i'=temp[_n+`i'] in `t'
				}
			}
			cap drop test*
			forvalues i=1(1)12 {
				quietly: gen test`i' = F`i'.sni-sni`i'
			}
		}
		}
		sum test*
	}
	}
*/
	
	
	*------------------------------------------------------
	* Forecasting monthly SNI
	*------------------------------------------------------

	* Estimating prefered model: ARIMA(1,0,1).  
	* 
	*   NOTE: This can also be motivated by referring to Donald W. K. Andrews and Werner Ploberger, "Testing for Serial Correlation Against an ARMA(1,1) Process," Journal of the American Statistical Association, Vol. 91, No. 435 (Sep., 1996), pp. 1331-1342

	cap drop time
	generate time =_n
	tsset time
	
	arima sni if year>=1991, arima(1,0,1)

	forvalues i=1(1)12{
		cap drop sni`i'
		generate sni`i'=.
	}

	sum time	
	local T=`r(max)'
	forvalues t=1(1)`T'{
	
		cap drop temp
		predict  temp, y dynamic(`t')
		
		forvalues i=1(1)12{
			replace sni`i'=temp[`t'+`i'] in `t'
		}
	}
	drop temp
	
	tsline sni sni1 sni12

	
	* Inspect forecast errors
	
	cap drop test*
	forvalues i=1(1)12 {
		quietly: gen test`i' = F`i'.sni-sni`i'
	}
	sum  test*
	drop test*


	
	*------------------------------------------------------
	* Create forecast errors = news shocks monthly SNI
   /*------------------------------------------------------	
   
	Note: The public formula is
	
			D_t = (G_t - X_t)/N_t, where
			
			D_t is the PF dividend in (fiscal) year t
			G_t is the gross distribution before expenses X_t
			N_t is the (estimated) number of eligible applicants 
			
			G_{t,m} = 0.5 * 21% * [\sum_{m=7}^{6} SNI^monthly_{t,m} + \sum_{a=4}^{1} SNI^{annual}_{t-a}],
			
			where m is the month in the current fiscal year t (i.e. starting from July of the previous 
			calendar year up to June of the current calendar year).	

			Therefore, the monthly news shock is the change in the forecast of 
			the monthly SNI times (0.5 x 21% / N_t).
	*/
	
	tsset Month

	
	* transform back to nominal forecasts to apply public dividend formula

	forvalues i=1(1)12{
		replace sni`i' = sni`i'*cpiu/100
	}
	replace sni = SNI_nominal


	* forecast current year's SNI
	
	cap drop SNIpredicted 
	generate SNIpredicted =.
	lab var  SNIpredicted "current fiscal year's predicted SNI"

	cap drop accumulated
	generate accumulated = 0

	local month = 6
	cap drop predicted
	egen predicted = rowtotal(sni1-sni12) if month==`month'
	replace accumulated = 0 if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 7
	cap drop predicted
	egen predicted = rowtotal(sni1-sni11) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 8
	cap drop predicted
	egen predicted = rowtotal(sni1-sni10) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 9
	cap drop predicted
	egen predicted = rowtotal(sni1-sni9) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 10
	cap drop predicted
	egen predicted = rowtotal(sni1-sni8) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 11
	cap drop predicted
	egen predicted = rowtotal(sni1-sni7) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 12
	cap drop predicted
	egen predicted = rowtotal(sni1-sni6) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'
	
	local month = 1
	cap drop predicted
	egen predicted = rowtotal(sni1-sni5) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 2
	cap drop predicted
	egen predicted = rowtotal(sni1-sni4) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 3
	cap drop predicted
	egen predicted = rowtotal(sni1-sni3) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 4
	cap drop predicted
	egen predicted = rowtotal(sni1-sni2) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 5
	cap drop predicted
	egen predicted = rowtotal(sni1-sni1) if month==`month'
	replace accumulated = accumulated[_n-1] + sni  if month==`month'
	replace SNIpredicted = accumulated + predicted if month==`month'

	local month = 6
	cap drop SNIrealized
	generate SNIrealized = accumulated[_n-1] + sni  if month==`month'
	replace  SNIrealized = 1030.5 if year==1991 & month==`month'

	drop predicted accumulated
	drop sni1-sni12


	* add number of eligible applicants, N_t
	
	preserve
	
		import excel "$homedir/data/APF_Data.xlsx", sheet("Eligible Applic + Div Adj.") cellrange(A2:I39) clear

		keep A G H I

		rename A FY // fiscal year
		rename G PFD
		rename H EligibleApplicants
		rename I DistributionAdjustment

		tempfile Applicants
		sort FY
		save `Applicants'
		
	restore

	generate FY = year   if month<=6 // generate fiscal year
	replace  FY = year+1 if month>=7
	bysort FY: generate FYmonth = _n
	order    FY FYmonth, before(year)

	replace FY = 1990 if year==1990 & FY==1991

	merge m:1 FY using `Applicants'
	drop if _merge==2
	drop _merge
	

	* Lag shock by one month since the numbers are end-of-month.
	*
	*   -> For instance, the August report is of August 30. Therefore, 
	*      the change between August and July reports is the relevant 
	*      shock for consumption response during September

	tsset Month

	foreach var in SNIpredicted SNIrealized sni SNI_nominal PFD EligibleApplicants DistributionAdjustment {
	
		cap drop `var'_orig
		generate `var'_orig = `var'
		
		cap drop temp
		generate temp = `var'
		
		replace `var' = L.temp
	}
	drop temp


	* calculate forecast error, G_t - G_{t-1}
	
	tsset Month
	cap drop DeltaG
	generate DeltaG = D.SNIpredicted if month!=7
	replace  DeltaG = SNIrealized - SNIpredicted[_n-1] if month==7 // NOTE: This makes the forecast errors predictable, since June also contains the adjustments to the SNI, but it also contains useful information

	
	* news shock = (0.5 * 21% / N_t) * (G_t - G_{t-1})
	
	cap drop news_sni
	generate news_sni = 0.5 * (21/100) / EligibleApplicants * DeltaG * 1000000
	lab var  news_sni "dividend news shock"
	
	lab var  sni "fund's net income"
	
	twoway connected news_sni Month, lc(blue) lw(medium) m(i) mlw(thin) mc(blue)  connect(stepstair) yaxis(1) || ///
	       connected sni      Month, lc(gray) lw(vthin)  m(i) mlw(thin) mc(black) connect(stepstair) yaxis(2) || ///
		   , xtitle("") graphregion(color(white)) legend(position(12) ring(0))
	
	corr news_sni D.sni sni L.sni

	
	* test whether forecast is biased
/*	
	cap drop test
	generate test = news_sni
	replace  test = test/cpiu*100
	reg test L(1/5).test, robust
	
	sum test, detail
	
	cap drop DeltaGres
	predict  DeltaGres, res
	sum DeltaGres DeltaG, d
	
	corr L(0/5).test DeltaGres
	
	
	* construct "unbiased" news shock
	
	tsset Month
	reg news_sni L(1/10).news_sni
	
	ac  news_sni if year>=1991, level(99)  // correlogram -> suggest q=1 or 6
	pac news_sni if year>=1991, level(99) // partial correlogram -> suggests p=2 or 5
	* => arima(1,0,1) seems resonable

	arima news_sni, arima(5,0,6)
	cap drop news_sni_res
	predict  news_sni_res, res
	lab var  news_sni_res "residualized news shock based on fund's SNI"
	
	reg news_sni_res L(1/10).news_sni_res
	
	tsline news_sni_res news_sni
*/
	
	
	*------------------------------------------------------
	* Construct expected dividend for figure
	*------------------------------------------------------
	
	* add previous 4 annual SNI for 5-year average

	preserve

		import excel "$homedir/data/APF_Data.xlsx", sheet("Monthly SNI") cellrange(A2:F310) clear

		keep A B F
		
		rename A year
		rename B month
		rename F sni_annual

		destring year,  replace
		destring month, replace

		generate FY = year   if month<=6
		replace  FY = year+1 if month>=7
		
		replace  FY = year if year<=1989
			
		drop if sni_annual==.

		tsset FY
		forvalues i=0(1)4 {
			generate SNI_`i' = L`i'.sni_annual
		}

		keep if FY>=1981
		sort FY
		tempfile SNIannual
		keep FY SNI_*
		save `SNIannual', replace

	restore

	merge m:1 FY using `SNIannual'
	drop if _merge==2
	drop _merge


	* sum SNI in previous 4 fiscal years

	egen SNIpast = rowtotal(SNI_1-SNI_4)
	
	
	* construct expected dividend
		
	cap drop div_exp
	generate div_exp = ( (SNIpast + SNIpredicted)/2 * 0.21 + DistributionAdjustment ) *1000000 / EligibleApplicants
	
	replace  div_exp = PFD_orig[_n-1] if month==7
	replace  div_exp = PFD_orig[_n-2] if month==8
	replace  div_exp = PFD_orig[_n-3] if month==9
	replace  div_exp = PFD_orig[_n-4] if month==10
	lab var  div_exp "expected PFD (market-based)"
	
	cap drop PFDactual
	generate PFDactual = div_exp // PFD
	replace  PFDactual = . if month!=10
	lab var  PFDactual "actual Permanent Fund Dividend (PFD)"
	
	replace  div_exp =. if month==10 
	
	cap drop PFDdash
	generate PFDdash = PFDactual
	forvalues i=1(1)12 {
		replace  PFDdash = PFDdash[_n+1] if PFDdash==.
	}
	lab var  PFDdash "actual Permanent Fund Dividend (PFD)"
	
	twoway connect div_exp   Month if Month>=m(1991m8) & Month<=m(2014m10), lc(black) lw(thin) m(i) mlw(thin) mc(black) connect(stepstair) cmissing(no)  || ///
		   connect PFDdash   Month if Month>=m(1991m8) & Month<=m(2014m10), lc(blue)  lw(vthin)  m(i) mlw(thin) mc(blue)  connect(stepstair) lp(shortdash) || ///	  		   
		   scatter PFDactual Month if Month>=m(1991m8) & Month<=m(2014m10), lw(vthin) m(o) msize(small) mlw(thin) mc(blue) || ///	  		   
		   , xtitle("") ytitle("dollars") yscale(r(0)) ylabel(0(500)2000) graphregion(color(white)) legend( position(7) ring(0) order(2 1) rows(2) region(lwidth(none)) ) 
	
		local file "$homedir/results/figures/dividend_vs_expectation_sni"
			graph export "`file'.eps", replace
			cap rm       "`file'.pdf"
			!epstopdf    "`file'.eps"
			graph export "`file'.tif", replace

		   

		*------------------------------------------------------
		* Assess forecast error
		*------------------------------------------------------
		   
		cap drop  ForecastError
		generate  ForecastError = PFDdash - div_exp
		summarize ForecastError, detail
		summarize ForecastError if Month<m(2008m1) | Month>m(2008m5), detail // exclude error due to Alaska Resource Rebate (to make forecast error comparable to the narrative-based forecast error) 

log close PFD_03
